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An  Introduction  to  the  Finite  Element  Method 


The  finite  element  method  views  a  structure  as  an  assemblage  of 
structural  elements  interconnected  at  a  finite  number  of  node  points. 
Consider  the  plane  elasticity  problem  in  Figure  1. 


We  wish  to  determine  both  the  deformation  and  stress  fields  for  the  region 
loaded  and  constrained  as  shown. 

The  essence  of  the  finite  element  method  is  as  follows: 

(1)  subdivide  the  continuum  into  a  number  of  small  simple 
u  elements,  as  shown  in  Figure  1 

(2)  assume  a  form  of  the  displacement  function  for  each  element 

type  (e.g. ,  linear,  quadratic,  etc.) 

(3)  derive  the  stiffness  relationship  for  each  element  (i.  e. ,  find 

Ke  in  Keu  =  f) 

(4)  assemble  the  element  stiffness  matrices  Ke  into  the  global 

stiffness  matrix  K  (this  is  easy  since  stiffness  can  be  added 
algebraically) 

(5)  solve  global  matrix  equation  Ku  =  F  for  the  displacement 

vector  u  (i.e. ,  the  components  of  u  are  the  displacements 
at  the  grid  points) 

(6)  knowing  the  displacements  everywhere,  compute  stresses. 
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The  above  will  be  illustrated  by  a  simple  example.  We  wish  to  derive 
a  triangular  membrane  element  which  is  suitable  for  plane  elasticity 
problems  (such  as  in  Figure  1). 

Consider  the  element  and  coordinate  system  shown  in  Figure  2. 


The  vertices  of  the  triangle,  called  grid  points  or  nodes,  can  displace 
in  either  the  x-  or  y- direction-  Thus,  the  element  shown  has  a  total 
of  six  degrees  of  freedom  (two  per  node). 

We  assume  displacement  functions  of  the  form 


u  =  a1  +  a2  x  +  ag  y 
v  =  a4*a5x  +  a6y 


where  u  and  v  are  the  displacement  components  in  the  x  and  y 
directions,  respectively.  Define  the  vectors  a  and  u  by 


where  u.  and  v.  are  the  u  and  v  components  at  point  i  . 
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Writing  equation  (1)  at  each  node  yields 

u  =  y  a 
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The  inverse  of  y  exists  and  is  given  by 
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The  strain  components  of  interest  will  be  grouped  into  a  strain 
vector  c  : 


For  this  example, 


e  ^ 

du/3x  j 

XX 

j 

dv/3y  | 

yy 

7xy 

5u/dy  +  dv/dx  J 

a2 

a6 

a3  +  a; 


(5) 


(6) 


(7) 
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(8) 


Relating  the  strain  vector  c  to  the  vector  a  defines  B: 

(  =  B  a 

where,  in  this  example, 

B  = 

Thus, 

Corresponding  to  the  strain  vector  is  the  stress  vector 


For  linear,  elastic  materials,  the  stress  and  strain  vectors  are  related 
by  Hooke's  Law: 

a  =  Die  (12) 

where  D  is  a  matrix  of  material  constants.  For  example,  for  plane 
stress  isotropy^  D  is  given  by 

‘1  v  0 

D  =  - 2  v  1  0  (13) 

l~V  0  0  (\-v)/2 

where  E  and  u  are  the  Young's  modulus  and  Poisson's  ratio  for  the 
material. 

Thus,  from  Equations  (10)  and  (12), 

a  =  D  B  a  u  (14) 


* 
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To  summarize,  Equation  (14)  gives  a  formula  for  the  stress  components 
for  an  element  given  the  displacements  u  of  the  node  points.  The 
matrices  D  ,  B,  and  a  depend  only  on  material  properties  and  geometry. 

The  final  step  for  the  element  is  to  compute  the  element  stiffness 
matrix  k.  The  general  result  will  be 

k  =  J  aT  BT  D  B  q  dV  (15) 

V 

where  the  integration  is  performed  over  the  volume  of  the  element. 

To  derive  (15),  let  P  be  the  vector  of  forces  at  the  nodes.  Applying 
an  arbitrary  virtual  displacement  6u  at  the  nodes,  (10)  yields 


6  (_  =  B  a  6  u  (16) 

By  the  principle  of  virtual  work,  the  work  done  by  P  at  the  nodes  must 
balance  the  internal  dissipation  of  energy;  thus, 

6uT  P  =  ;  6eT  ?dV  (17) 

V 


Substituting  (16)  and  (14)  into  (17)  yields 


6uT  P  =  '  6uT  aT  BT  D  B  o  u  d  V 


=  tiuT  (  ’ttTBTDBttdV)u 


(18) 


For  (18)  to  hold  for  an  arbitrary  variation  of  displacement  6u  , 

P  =  (  J  aT  BT  D  BodV)u  (19) 

By  definition,  the  parenthetical  expression  in  (19)  represents  the  stiffness 
matrix  for  the  element;  hence,  (15)  is  proved. 

For  the  triangular  membrane  element,  all  matrices  in  the  integrand 
of  (15)  are  constant,  so  that  k  becomes 

k  =  |  x2  y3  t  (gTBTDBg)  (20) 
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where  t  is  the  membrane  thickness.  This  matrix  is  shown  on  page  3-7. 

Note  that  equations  (8),  (10),  (12),  (14),  and  (15)  are  general,  so 
that  they  apply  to  other  finite  elements  as  well  as  to  the  triangular 
membrane. 

For  a  problem  with  many  elements,  the  stiffness  matrix  k  for  each 
element  can  be  computed  and  assembled  into  a  global  stiffness  matrix 
K  for  the  entire  structure.  Then,  one  must  solve  the  system 

Ku  =  F  (21) 


where  u  is  now  the  vector  of  nodal  displacements  for  the  entire  structure, 
and  F  is  the  vector  of  loads  applied  at  the  nodes.  Equation  (21)  is 
basic  to  the  finite  element  method. 

It  is  apparent  that  (21)  is  merely  a  generalization  of  the  familiar 
expression  for  a  one-dimensional  spring: 

kx  a  f  (22) 

Thus,  it  is  not  too  surprising  that,  for  time-dependent  problems,  the 
analogy  still  holds  and  one  obtains  the  matrix  equation 

Mu+Bu+Ku  =  F(t)  (23) 


Here,  and  B  are  mass  and  damping  matrices,  respectively,  and  the 
dots  over  u  indicate  differentiation  with  respect  to  the  independent 
variable  t  (time). 

In  essence,  the  approach  outlined  above  approximates  the  displacement 
solution  by  a  set  of  piecewise  polynomials.  It  can  thus  be  shown  to  be  a 
variant  of  the  Rayleigh-Ritz  method. 
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Stiffness  Matrix  for  Plane  Stress  Triangular  Membrane 
Finite  Element 


DEMONSTRATION  PROBLEMS 


See  "NASTRAN  Sample  Problem  Computer  Output"  by  G.C.  Everstine  and 
M.M.  Hurwitz  (DTNSRDC/CMLD-81-04)  for  problem  solutions. 


1.  Cantilever  Beam  with  Point  Load 


For  steel,  E  =  30x10®  psi ,  v  =  0.3,  p  =  7.324xl0~\b-sec^/in^ 
Find  displacements,  stresses,  and  reactions. 


1A.  Arch  Under  Static  Pressure 


Subcase  1 
Subcase  2 
Subcase  3 


Material:  Steel 

Pressure  Load 
Gravity  Load 

Solve  a  ring-stiffened  cylinder  problem  by  changing 
B.C.  at  0  =  t  40°  to  planes  of  symmetry. 


IB.  Ring-Stiffened  Cylinder  with  Pressure  Load 
(Modeled  using  conical  shell  elements) 

Find  displacements  for  same  cylinder  as  in  problem  1A,  Subcase  3. 


Find  displacements  by  modeling  only  one-half  of  structure. 


Linear  Steady-State  Heat  Conduction 


For  the  arch  structure  of  problem  1A,  find  the  steady-state 
temperature  distribution  due  to  a  uniform  applied  heat  flux  at  9  =  -  4CT  of 
0.0025  BTU/sec-in^,  with  the  edges  at  2  =  0  and  Z  =  96  immersed  in  an  ice 
bath  (32°F).  The  thermal  conductivity  of  steel  is  7.175x10*4  BT'J/sec-in-°F. 


IE. 


2-D  Poisson  Equation  (Torsion  of  Triangular  Prism) 


Compute  the  maximum  shear  stress  in 
a  twisted  bar  whose  cross-section  is 
an  equilateral  triangle  of  altitude 
a  =  0.09  m.  The  bar  has  a  shear 
modulus  G  =  80  GPa  and  is  subjected 
to  an  angle  of  twist  per  unit  length 
6  ■  0.04  rad/m.  (The  exact  solution 
for  the  maximum  shear  stress  is  Gea/2.) 


Find  natural  frequencies  and  mode  shapes  for  beam  of  problem  1. 
(INV,  GIV,  FEER,  and  eigenvalue  APPEND) 
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4.  Pi fferential  Stiffness 


5.  Buckling 

For  the  beam  of  problem  4,  at  what  value  of  F  and  M  (assuming  F 
and  M  are  numerically  equal)  will  buckling  occur? 

6.  Piecewise  Linear  Analysis 

Find  the  displacements  and  stresses  for  the  following  frame: 


y 


7.  Complex  Eigenvalue  Analysis  (Direct  Method) 

Find  the  damped  natural  frequencies  for  the  beam  of  problem  1 
with  a  dashpot  (c  =  0.25  lb-sec/in)  connected  in  the  y-direction 
between  point  13  (at  x  =  60)  and  ground. 

8.  Frequency  Response  (Direct  Method) 


1  75"  ■■  1  — ^F(t)  =  Fq  cos  ut 

F„  =  1  lb  ,  u>  =  2uf 

i 

Find  steady-state  displacement  response  at  f  =  3  Hz  and  7  Hz  for 
beam  of  problem  1. 
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9.  Transient  Response  (Direct  Method) 

Find  time-dependent  response  for  beam  of  problem  #8  with  F(t) 
given  below: 


Problem  9A  restarts  from  t= 100  to  illustrate  TRD  CONTINUE. 

10.  Complex  Eigenvalues  (Modal  Method) 

Same  as  problem  #7,  except  use  modal  approach. 

> 1 . &  1 1  A.  Frequency  Response  (Modal  Method) 

Same  as  problem  #8,  except  use  modal  approach.  (INV  and  GIV) 

1 2 .  Transient  Response  (Modal  Metnod) 

Same  as  problem  #9,  except  use  modal  approach. 

13.  Normal  Modes  with  Differential  Stiffness 

For  the  beam  of  problem  #1,  find  the  flexural  natural  frequencies 
for  the  beam  spinning  about  the  y  axis  at  1.5  Hz. 
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WORKSHOP  PROBLEMS 


material:  steel  (E 


30x1 06  psi ,  v  =  0.3,  p  = 


7 . 324x1 0~4  lb-sec2/in4) 


This  structure  will  be  used  for  all  types  of  analysis. 

Suggested  F.E.  mesh:  6x4  mesh  of  plate  elements  covering  the  30"xl2" 

region  with  beam  stiffeners 

(Note  existence  of  plane  of  symmetry.) 


1.  Static  Stress  Analysis.  Determine  stresses  and  displacements  for 

a.  uniform  unit  pressure  load  on  plate  in  -z  direction 

b.  gravity  load  in  -z  direction 

c.  sum  of  (a)  and  (b) 


la. 


Plotting. 

a.  Plot  undeformed  structure. 

b.  Plot  static  deformation  for  pressure  load  of  problem  #1. 

c.  Make  matrix  topology  plot  of  constrained  stiffness  matrix  (KLL) 
for  problem  -1 . 


lb. 


DMAP.  For  problem  #1,  starting  with  the  constrained  stiffness  matrix 
( KLL)  and  load  vectors  (PL),  write  a  DMAP  ALTER  to  R.F.  1  to: 

a.  compute  displacement  vectors  (to  be  called  UTEST) 

b.  compute  the  residual  (RESID  =  PL  -  KLL  x  UTEST) 

c.  compute  error  ratio 


ERR1 


RESID  •  UTEST 
PL  •  UTEST 
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d.  print  UTEST,  RES I D ,  and  ERR1  and  check  with  NASTRAN's  values. 

e.  do  problem  lc  both  before  and  after  grid  point  resequencing 
to  reduce  matrix  wavefront. 

2.  Inertia  Relief.  The  plate  is  free-free  (for  this  problem  only)  rather 
than  cantilevered.  Determine  stresses  resulting  from  a  100  lb.  normal 
point  force  at  the  center  of  plate. 

3.  Normal  Modes.  Find  the  lowest  3  natural  frequencies  and  modes. 

4.  Differential  Stiffness.  The  load  is  a  uniform  downward  (-z)  pressure  of 
2  psi  plus  a  20,000  lb  uniform  line  load  at  the  free  end  in  the 

-x  direction.  Determine  the  displacements. 

5.  Buckling.  Find  the  factor  by  which  the  load  of  problem  #4  should  be 
scaled  to  induce  buckling. 

6.  Piecewise  Linear  Analysis.  Determine  the  stresses  and  displacements  for 
a  uniform  downward  10  psi  pressure  load  on  the  plate.  Assume  that  steel 
has  the  symmetric  stress-strain  relation: 


7.  Complex  Eigenvalues.  Add  uniform  viscous  dampers  in  the  z-direction 
which  are  attached  between  the  plate  and  ground  along  y=3  and  y=9.  The 
total  damping  along  each  line  is  2.4  Ib-sec/in.  Find  the  damped 
frequencies  and  modes. 

8.  Frequency  Response.  Compute  frequency  response  at  the  free  end  to 
sinusoidal  load  at  the  free  end  of  each  stiffener:  F  =  10  cos  2Tift, 
where  f  =  30,  32,  34,  ....  44. 

8a.  Random  Response.  Determine  the  power  spectral  density  (PSD)  function 
and  rms  value  of  the  z-di spl acement  at  the  center  of  the  free  end 
(x=30,  y =6 )  for  the  input  normal  force  FSD  (relative  to  the  loading  of 
problem  *8)  at  the  free  end  of  each  stiffener  given  on  the  following 
page. 
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PSO  ,  lb2 -sec 


Plotting.  Make  XY  printer  plots  vs.  frequency  of  the  responses 
computed  in  problems  #8  and  8a. 

Transient  Response.  Determine  the  first  75  msec  of  response  at  the 
free  end  due  to  the  following  normal  load  at  the  end  of  each 
stiffener: 


Plotting.  Make  XY  plot  (both  on  paper  and  on  plotter)  for  transient 
z-displacement  response  at  center  of  free  end. 


BANDIT  EXAMPLE 


BEFORE  '  AFTER 

(MATRIX  BANDWIDTH  =  63)  (MATRIX  BANDWIDTH  =  17) 


TIME  AND  CORE  ESTIMATION 


The  general  approach  to  estimating  NASTRAN  CPU  run  times  is  to 
identify  those  functional  modules  which  will  be  most  time-consuming, 
estimate  the  CPU  time  for  those  modules,  and  then  assume  that  the  total 
CPU  run  time  equals  about  140%  of  the  sum  of  the  estimated  times. 

Every  rigid  format  generally  requires  stiffness  matrix  and  mass  matrix 
generation  and  at  least  one  matrix  decomposition.  The  methods  for  computing 
these  times  will  be  given  first.  Other  details  for  each  rigid  format  will 
then  be  given. 

Core  requirements  are  usually  thought  of  as  the  core  required  to  have 
a  "no  spill"  condition  for  decomposition.  These  requirements  will  be 
discussed  in  the  Matrix  Decomposition  section. 

Matrix  Generation 

The  functional  modules  for  stiffness  and  mass  matrix  generation  are 
SMA1  and  SMA2  for  Level  15  (and  below)  and  EMG  for  Level  16  and  above.  The 
CPU  times,  for  various  machines,  for  the  generation  of  one  element  stiffness 
matrix,  for  various  elements,  are  given  in  Table  7-1  for  NASTRAN  Levels  15, 
16,  and  17.  The  time  required  to  compute  an  element  lumped  mass  matrix  is 
negligible,  while  the  time  to  compute  an  element  consistent  mass  matrix  is 
approximately  the  same  as  the  time  required  for  the  element  stiffness  matrix. 

Matrix  Decomposition 

The  time  required  to  decompose  a  real ,  symmetric  matrix  with  no  spill  is 
T  =  ^M(nb2  -  |b3)  (Level  15  and  below)  (7-1) 

T  =  j  MnC^s  (Level  16  and  above)  (7-2) 

where 

T  =  CPU  time 

n  =  order  of  the  matrix 

b  =  matrix  semi -bandwidth 

C _  =  matrix  rms  wavefront 

rms 

M  =  machine  time  constant  (see  Table  7-3) 


Table  7-1. 

element  matrix  generation  times  for  nastran 


No .  of 

GDC  6400 

Element 

Gauss  Pts . 

Time  (Sec 

)•'  Remarks 

BAR 

0.27 

BEAM 

0.22 

MSC/NASTRAN 

CONEAX 

0.43 

per  harmonic 

HEXA1 

2.0 

HEXA2 

4.0 

HEXA  (8  nodes) 

0.80 

MSC/NASTRAN 

HEXA  (20  nodes) 

5.60 

MSC/NASTRAN 

HEX20 

6.40 

MSC/NASTRAN 

IHEX1 

2 

1.0 

3 

2.5 

4 

5.5 

2 

0.2 

heat  conduction 

IHEX2 

2 

5.0 

3 

12.5 

3 

4.0 

single  precision  (DTNSRDC) 

4 

27.5 

I  HEX  3 

3 

31.0 

4 

68.0 

IS2D8 

3 

2.2 

single  precision  (DTNSRDC) 

IS3D8 

2 

2.4 

Sperry/NASTRAN 

IS3D20 

2 

13.4 

Sperry/NASTRAN 

3 

39.2 

Sperry/NASTRAN 

PENTA  (6  nodes) 

0.33 

MSC/NASTRAN 

PENTA  (15  nodes) 

1.80 

MSC/NASTRAN 

QDMEM 

0.60 

QU ADI, QUAD 2 

2.02 

QUAD  4 

0.30 

MSC/NASTRAN 

ROD 

0.06 

SHEAR 

0.30 

TETRA 

0.40 

TRAP AX 

2.20 

TRAPRG 

1.05 

TRIAAX 

0.55 

TRIA1 , TRIA2 

1.52 

TRIA3 

0.14 

MSC/NASTRAN 

TRIARG 

0.63 

TRIAX6 

1.05 

MSC/NASTRAN 

TRIM6 

0.28 

TRMEM 

0.15 

TUBE 

0.06 

WEDGE 

1.20 

*See  Table  7-2  for  conversion  factors  for  other  computers. 
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Table  7-2. 

CONVERSION  FACTORS  FOR  OTHER  COMPUTERS 


To  Convert  the  CDC  6400 

Times  in  Table  7-1  to 

Multiply  by 

CDC  6500  and  Cyber  73 

1.0 

CDC  6600  and  Cyber  74 

0.33 

CDC  7600  and  Cyber  76 

0.05 

Cyber  172 

0.65 

173,  174 

0.45 

175 

0.12 

176 

0.05 

IBM  360/50 

2.0 

65 

0.60 

75 

85 

0.18 

91,  95 

0.15 

370/155 

0.75 

158 

165 

0.18 

168 

195 

Univac  1108 

0.33 
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The  timing  equations  (7-1)  and  (7-2)  are  just  the  expected  dominant 

terms  of  the  full  timing  equations,  which  may  be  found  in  the  Theoretical 
Manual.  For  Equation  (7-1),  the  number  of  active  columns  outside 

the  band  (in  the  Level  15  sense)  is  assumed  to  be  small,  as  would  normally 
be  the  case  if  a  bandwidth  resequencing  program  is  used. 

The  core  required  to  perform  a  real,  symmetric  decomposition  with  no 
spill  is 

N  =  j  b(b+5)p  +  PROG  (Level  15  and  below)  (7-3) 

N  =  i  p  +  PROG  (Level  16  and  above)  (7-4) 

c  mdx 


where 

N  =  the  number  of  decimal  words  required 
b  =  matrix  semi -bandwidth 
p  =  1  for  CDC,  2  for  IBM  and  UNI VAC 

PROG  (given  in  decimal  words)  =  20000  for  CDC,  30000  for  IBM, 

. .  25000  for  UN  I  VAC 

C  =  matrix  maximum  wavefront 
max 

The  first  term  on  the  right-hand  sides  of  Equations  (7-3)  and  (7-4)  may  be 
termed  the  working  storage  W. 

The  Level  15  formula.  Equation  (7-3),  assumes  a  small  number  of  active 
columns. 

The  bandwidth,  b,  and  wavefront  terms,  C  and  C  ,  required  in 

rms  max 

equations  (7-1)  through  (7-4)  may  be  computed,  using  BANDIT,  as  follows: 

Multiply  the  required  term  (b,  C  ,  ,  or  C  )  obtained  from  BANDIT  by 

max  rms 

6  _  (MPC  +  SPC  +  PS) 


where 


SPC  =  number  of  degrees  of  freedom  SPC’d 

MPC  =  number  of  dependent  degrees  of  freedom  MPC'd 


PS  =  number  of  degrees  of  freedom  constrained  on  GRID  cards 

G  =  number  of  grid  points  in  the  problem 

(This  formula  represents  the  average  number  of  free  degrees  of  freedom  per 

grid  point.)  If  a  significant  number  of  degrees  of  freedom  are  OMIT'd, 

then  it  may  be  assumed  that  b=n,  C  =n,  and  C  =  n//3.  (n=matrix  order) 

J  max  rms 

For  real ,  unsymmetric  decomposition  with  no  spill,  for  both  Levels  15 

and  16,  the  time  and  working  storage  required  are  each  p  times  the  amount 

required  for  real,  symmetric  decomposition  on  Level  15,  where  p=4  for  CDC, 

2  for  IBM  and  UNIVAC. 

For  Level  16  complex,  symmetric  decomposition  with  no  spill,  the  time 
is  4  times  as  long  as  Level  16  real,  syrmetric  decomposition,  and  the 
working  storage  must  be  twice  that  required  for  Level  16  real,  symmetric. 

For  Level  15  complex  decomposition,  both  symmetric  and  unsymmetric, 
and  for  Level  16  complex ,  unsymmetric  decomposition,  the  time  is  16-20  times 
longer  than  Level  15  real,  symmetric  decomposition.  The  working  storage 
required  is  p  times  the  amount  required  for  Level  15  real,  symmetric 
decomposition,  where  p=8  for  C0C,4  for  IBM  and  UNIVAC. 

Forward-Backward  Substitution  (FBS) 

Forward -backward  substitution  (FBS)  is  the  second  (and  final)  step 
required  in  the  solution  of  a  set  of  simultaneous  linear  equations.  (Matrix 
decomposition  is  the  first  step.)  Normally,  FBS  times  are  relatively  small. 
However,  it  is  not  uncommon  for  FBS  times  to  be  significant.  The  time 
required  for  FBS  is: 

T  =  2bnrM  +  [-£J-]  2nblu  (Level  15)  (7-5) 

T  =  2nCavg(rM  +  [-^-]  P$)  (Level  16)  (7-6) 

where 

r  =  number  of  right-hand  side  vectors 

Cflvg  =  average  wavefront  (determined  from  BANDIT) 

W  =  working  storage 


Once  again,  the  terms  in  Equations  (7-5)  and  (7-6)  are  the  dominant  ones  in 
the  full  timing  equations.  Brackets  mean  "next  larger  integer". 


Matrix  Multiply  anc!  Add  (ITPVAD) 


A  B  =  D 
mn  np  rnp 


(7-7) 


T  =  mnpA  (p  M  +  [ nj)  *  3  P$  )  (Level  16)  (7-8) 

where  is  the  density  of  the  A  matrix  and  the  timing  constants  are 

given  in  Table  7-3. 


Multipoint  Constraints  (MPC) 


In  multipoint  constraint  elimination,  three  multiply-add  (MPYAD) 
operations  are  performed.  If  we  use  the  notation  of  Equation  (7-7)f 
then  the  matrix  orders  of  the  three  multiplications  are: 


1. 

2. 

3. 


m=d,  n=d,  p=i, 
m=i,n=d,  p=i, 
m=i,n=d,p=i, 


P  = 
P  = 
P  = 


density  of  stiffness  matrix  (BANDIT) 

density  of  stiffness  matrix  (BANDIT) 

density  of  multipoint  constraint 
transformation  matrix 


where: 


d  =  number  of  MPC  equations 

i  =  (total  number  of  degrees  of  freedom  in  the  problem)  -  d 


The  density  of  the  transformation  matrix  may  be  estimated  as  the  ratio  of  the 
number  of  degrees-of-freedom  in  MPC  equations  to  the  total  number  of  aegrees- 
of- freedom  in  the  problem. 

Guyan  Reduction 

OMIT  partitioning  and  Guyan  Reduction  may  be  time-consuming  processes. 
One  real,  symmetric  matrix  decomposition  is  required,  the  order  of  the 
matrix  being  the  number  of  degrees  of  freedom  in  the  o-set,  i.e.,  degrees 
of  freedom  omitted.  Also,  there  will  be  as  many  FBS's  as  there  are  degrees 
of  freedom  in  the  a-set,  i.e.,  degrees  of  freedom  not  MPC'd,  SPC'd,  or 
OfilT'd. 


For  some  of  the  more  commonly  used  rigid  formats,  additional  timing 
details  and  core  requirements  are  as  follows: 
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Rigid  Format  1 

In  static  analysis,  usually,  the  only  major  time-consumi ng  operations 
are  stiffness  matrix  generation,  real  symmetric  decomposition,  anu 
multipoint  constraints. 

Rigid  Format  3 

Time  for  one  eigenvalue  (Functional  Module  READ)  by  INV  =  one  decomposition 
+  8  F8S 

3 

Time  for  Givens  method  =  An  M,  where 

A  =  factor  of  5-10  depending  on  the  number  of  eigenvectors  requested, 

5  for  none,  10  for  all 
n  =  matrix  order 

M  =  multiply  times  as  previously  given 
The  core  required  for  Givens  method  is  computed  as  though  the  matrix  were 
full,  i.e.,  bandwidth  =  maximum  wavefront  =  n. 

Rigid  Formats  4  and  5 

The  differential  stiffness  matrix  time  is  approximately  the  same  as 
the  linear  stiffness  matrix  time. 

Rigid  Format  8 

Time  for  one  requested  frequency  (Functional  Module  FRRD)  is  approxi¬ 
mately  the  time  for  a  matrix  decomposition.  In  Level  15,  the  decomposition 
is  always  complex,  unsymmetric.  In  Level  16,  the  decomposition  depends  on  the 
matrix,  and,  therefore,  may  be  real  or  complex,  symmetric  or  unsyimetric. 

The  inclusion  of  DMIG  cards  will  automatically  trigger  an  unsymmetric 
decomposition. 

Rigid  Format  9 

Time  for  transient  integration  (Functional  Module  TRD)  =  one 

decomposition  per  time  step  change  +  two  FBS  per  time  step  (with  r=l 
in  Equations  (7-5)  and  (7-6)) 

Inclusion  of  DMIG  cards  will  automatically  trigger  an  unsynmetric 
decomposition. 
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Rigid  Formats  11  and  12 

The  dominant  time  in  these  rigid  formats  will  come  from  the  eigenvalue 
analysis  required  to  compute  the  requested  mode  shapes,  not  from  the  frequency 
response  or  transient  response  analyses.  Therefore,  the  time  will  be 
approximately  the  same  as  in  Rigid  Format  3. 

If  the  user  adds  the  times  for  stiffness  and  mass  matrix  computation  to 
the  times  given  for  these  rigid  formats,  and  increases  the  rotal  time  by  about 
40%,  the  result  should  be  close  to  the  total  NASTRAN  CPU  time. 
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BASIC  PLOTTING 


Types 

1.  structural  - 

undeformed  or  deformed 

orthographic,  perspective,  or  stereoscopic  projection 

2.  X-Y 

3.  topological  displays  of  matrices 

4 .  contour 

Structural  Plotting  (Theory) 

The  model  is  defined  in  the  basic  coordinate  system  (XY2). 

The  plotter  coordinate  system  is  denoted  RST.  The  XYZ  system 
(fixed  with  respect  to  the  structure)  is  oriented  with  respect 
to  the  RST  system  (fixed  with  respect  to  the  plotter)  in  two 
steps : 

1.  overlay  XYZ  on  RST  in  some  order 

2.  rotate  XYZ  w.r.t.  RST  by  the  angles  y,  6,  a  (in 
that  order),  which  are  the  angles  of  rotation  about  the  T,  S,  and 
R  axes,  respectively. 

For  orthographic  plots,  NASTRAN  then  plots  the  projection  in 
the  ST-plane.  For  perspective  plot:,  NASTRAN  also  needs  to  know 
the  location  of  a  vantage  point  and  of  a  projection  plane  (plotter 
surface).  Stereoscopic  plots  consist  of  two  perspective  images, 
each  with  a  slightly  different  vantage  point,  which  are  viewed 
simultaneously  using  a  stereo  viewer. 

Deformed  plots  also  require  the  user  to  specify  the  scaling 
to  be  applied  to  the  deformation. 

Structural  Plotting  (Practice) 

1.  All  NASTRAN  plots  are  written  on  a  file  called  PLT1  or 
PLT2  (usually  PLT2 ) .  Therefore,  using  system  control  cards,  user 
must  request  computer  operator  to  mount  a  plot  tape  of  that  name. 
At  some  installations,  plot  files  can  be  written  to  disk. 

2 .  Insert  plot  control  cards  into  the  case  control  deck 
immediately  before  the  Begin  Bulk  card.  (See  next  section.) 

3.  If  necessary,  request  that  the  plot  tape  be  processed 
(i.e.,  plotted)  by  the  plotter. 
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Plot  Control  Cards  (Structural  Plots) 

1.  Example  1:  undeformed  plot 

PLOTID  =  JOHN  DOE,  NSRDC  1844 
OUTPUT ( PLOT ) 

PLOTTER  SC  MODEL  4020 

SET  1=4,8,10  THRU  27 ,  41 

ORTHOGRAPHIC  PROJECTION 

AXES  Z,  X,  Y 

VIEW  -40 .  ,  30 .  ,  0  . 

FIND  SCALE,  ORIGIN  1,  SET  1 
PLOT  SET  1 

2.  Example  2:  deformed  plot 

Insert  the  following  before  FIND  card: 

MAXIMUM  DEFORMATION  3. 

Replace  PLOT  card  with 

PLOT  STATIC  DEFORMATION,  SET  1 

Plot  Control  Cards  (XY  Plotting) 


1.  Example  1: 

PLOTID  =  JOHN  DOE,  NSRDC  1844 
OUTPUT (XYPLOT) 

PLOTTER  =  SC  4020 
parameter  cards  (optional) 

XYPLOT  DISP  RESPONSE  2,5/16  (Tl) 
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Perspective 


NASTRAN  Plotting  at  DTNSRDC 


The  CAI.COMP  plotter,  model  936,  is  used  to  process  NASTRAN  plots 
at  DTNSRDC.  Since  this  model  is  not  supported  directly  by  NASTRAN,  the 
"plotter"  specified  by  the  user  is  NASTPLT,  the  general  purpose  plotter 
package.  The  plotting  output  written  by  NAS  IRAN  on  tile  PLIJ  is  then 
interpreted  by  a  post-processor  called  PLTTRN9  36 ,  which  was  written  by 
J.M.  McKee.  The  output  from  PLTTRN936  is  written  on  a  tape  (assigned  by 
the  user)  and  processed  by  the  CALCOMP  936. 

Thus,  the  PLOTTER  card  in  the  plot  request  packet  of  the  NASTRAN 
Case  Control  Deck  is 

PLOTTER  NASTPLT,  M0DEL  T,  0 

The  additions  to  the  CDC  system  control  cards  are  as  follows: 

1.  Prior  to  the  NASTRAN  execution  card,  insert 

REQUEST , PLT2 ,  -PF . 

2.  After  the  NASTRAN  execution  card,  insert 

CATALOG, PLT2, . . . 

REWIND, PLT2. 

ATTACH ,  PLTTRN ,  PLTTRN9  36  ,  IDCA.MY ,  MR=1 . 

UNLOAD, TAPE 3 . 

REQUEST, TAFE3, HI, S, RING, VSN=XXXXXX. 

PLTTRN. 

3.  The  MT  parameter  on  the  job  card  must  be  increased,  if 
necessary,  by  unity. 

In  order  for  PLT2  to  be  written  on  disk  rather  than  on  tape,  PLT2 
should  be  listed  after  the  FILES  parameter  of  the  NASTRAN  card,  which 
precedes  the  ID  card.  Thus,  for  example, 

NASTRAN  C0NFIG=6 , FILES=(NPTP ,0PTP, PLT2) 

A  plot  request  (form  10462/26)  must  be  submitted  to  ADP  Control  for 
each  tape. 


i 
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NASTRAN  ELEMENTS 


Elements  are  defined  on  connection  cards  (e.g.,  CBAR,  CROD; , 
which  list  the  grid  points  to  which  the  elements  are  connected 
and  refer  to  property  cards  (e.g.,  PBAR,  PROD),  which  define 
geometrical  properties  (e.g.,  A,  I,  J,  t,  etc.)  and  in  turn 
refer  to  material  property  cards  (e.g.,  MAT 1 ) ,  which  define 
material  properties  (e.g.,  E,  G,  v,  p). 

The  elements  contained  in  level  17  are  summarized  and 
categorized  in  Table  9-1  (p.  9-8)  and  described  (by  category)  as  follows: 

One-dimensional  elements 
CBAR  (beam) 

based  on  simple  beam  theory 
assumes  uniform  properties  over  the  length 
assumes  shear  center  coincides  with  elastic  axis 
includes  extension,  torsion,  bending  in  two  planes,  shear 
ends  may  be  offset  from  defining  grid  points 
any  5  of  6  forces  at  each  end  may  be  set  equal  to  zero 
using  pin  flags 
6  DOF/point 
CROD  (rod) 

special  case  of  beam  with  only  axial  and  torsional 
properties 

no  offsets  or  pin  flags 
2  DOF/point  (u  ,  R  ) 

X  X 

CONROD  (rod) 

same  as  ROD  except  properties  are  included  on  connection 
card 

CTUBE  (tube) 

a  rod  of  circular  cross-section,  either  solid  or  hollow 
CVISC  (viscous  rod) 

extensional  and  torsional  viscous  damping  properties  rather 
than  stiffness  properties 
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Two-dimensional  elements 

CTRMEM  (triangular  membrane) 

constant  strain  triangle  (CST) 

2  DOF/point  (u^,  u  ) 

CQDMEM  (quadrilateral  membrane) 

two  pairs  of  overlapping  TRMEM’s 

CQDMEMl  (quadrilateral  membrane) 

linear  isoparametric  quadrilateral  membrane  element 
2  DOF/point  (u^,  u  ) 

inefficient  since  a  4x4  Gauss  quadrature  is  performed 
instead  of  2x2 

most  accurate  of  the  quadrilateral  membranes  in  !!  ASTRA! 
CQDMEM2  (quadrilateral  membrane) 
four  non-overlapping  TRMEM's 

2  DOF/point  (u  ,  u  ) 

x  y 

more  accurate  than  QDMEM 

CTRBSC  (basic  bending  triangle) 

3  DOF/point  (ny,  R^,  R  ) 
normal  displacement  u^  varies  as  (incomplete)  cubic  in 

x  and  y 

used  as  building  block  for  other  elements  rather  than 
as  stand-alone 

CTRPLT  (triangular  bending  plate) 

three  non-overlapping  basic  bending 
triangles  joined  at  centroid  (the 
Clough  bending  triangle) 

3  DOF/point  ( u z ,  R^ ,  R  ) 

CQDPLT  (quadrilateral  bending  plate) 

two  pairs  of  overlapping  basic  bending  triangles 

3  DOF/point  (u  ,  R  ,  R  ) 
z  x  y 
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CTRIAi  (membrane  and  bending  triangle) 

superposition  of  TRMEM  and  TRPLT 

5  DOF/noint  (u  ,  u  ,  u  ,  R  ,  R  ) 
x  y  z  x  y 

i  =  1:  nonhomogeneous  panel  (e.g.,  sandwich  or  honeycomb 
construction 
i  =  2:  homogeneous  panel 

CQUADi  (membrane  and  bending  quadrilateral) 

superposition  of  QDMEM  and  3DPLT 

5  DOF/point  (u  ,  u  ,  u  ,  R  ,  R  ) 

^  x  y  z  x  y 

1  =  1,2:  same  as  for  TRIAi 

CTRIM6  (triangular  membrane) 

linear  strain  triangle  (LST) 
six  nodes  (three  corner,  three  mid-side) 
thickness  can  vary  bilinearly  in  x  and  y 

2  DOF/point  (u^,  u^) 

most  accurate  of  the  NAS TRAN  membranes 

CTRPLTl  (triangular  bending  plate) 

six  nodes  (three  corner,  three  mid-side) 

3  DOF/point  (u  ,  R  ,  R  ) 

z  x  y 

normal  displacement  u^  varies  as  (incomplete)  quintic 
in  x  and  y 

thickness  can  vary  bilinearly  in  x  and  y 
CTRSHL  (triangular  shallow  shell) 

a  combination  of  TRIM6  and  TRPLT1 

shell  surface  (which  need  not  be  flat)  is  approximated 

quadratically  in  x  and  y 

5  DOF/point  (u  ,  u  ,  u  ,  R  ,  R  ) 
x  y  z  x  y 

CSHEAR  (quadrilateral  shear  panel) 

resists  the  action  of  tangential  forces  applied  to  its 

edges  but  not  the  action  of  normal  forces 

usually  used  in  combination  with  rods  or  beams 

1  DOF/point  (m-plane  along  diagonal) 


l'TWIST  (quadrilateral  twist  panel) 
bending  analog  of  shear  panel 

equivalent  for  bending  action  to  a  pair  of  parallel 
shear  panels 

1  DOF/point  (moment  having  axis  perpendicular  to  diu jonai 
and  in-plane) 


Three-dimensional  elements 


CTETRA  (tetrahedron) 

constant  strain  tetrahedron 
3-D  analog  of  TRMEM  (CST) 

4  vertices,  4  triangular  faces 
3  DOF/point  (u^,  u  ,  u^) 

CHEXAi  (hexahedron) 

8  vertices,  6  quadrilateral  faces 
superposition  of  5  non-overlapping  TETRA 
elements  (HEXA1)  or  10  overlapping 
TETRA  elements  (HEXA2) 

3  DOF/point  (u^,  u^ ,  uz) 


CWEDGE  (wedge) 

6  vertices,  3  quadrilateral  faces  and 
2  triangular  faces 
superposition  of  3  TETRA1 s 
3  DOF/point  (u^,  u^ ,  uz) 

CIHEXi  (isoparametric  hexahedron) 

8  vertices,  6  faces 

linear  (i  =  l)  ,  quadratic  (i  =  2)  ,  or  cubic 
(i=3)  shape  functions 
8  (i~l),  20  (i=2),  or  32  (i=3)  nodes 
isotropic  materials  only 
3  DOF/point  (u  ,  u  ,  u  ) 


i 

1 
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Axisymmetric  Elements 


CCONEAX  (conical  shell) 

for  axisymmetric  thin  shells 

includes  membrane,  bending,  and  transverse 

shear  effects  n 

loads  may  be  non-axisymmetric  since  \ 

motions  are  expanded  in  Fourier  . \ 

series  in  aximuthal  coordinate  _ 

5  DOF/point  (normal  rotation  is  excluded) 

CTORDRG  (doubly  curved  toroidal  ring) 
for  axisymmetric  thin  shells 

axisymmetric  loads  only  , 

includes  membrane  and  flexural  behavior 

membrane  displacement  function  is  complete  \ 

cubic  _ _ 

flexural  displacement  function  is  complete 
quintic 

5  DOF/point  (u  is  excluded) 

U 

CTRIARG  (triangular  ring)  z 

solid  of  revolution  element  for  thick- 

walled  axisymmetric  structures  P'S. 

loads  must  be  axisymmetric  l  / 

linear  displacement  function 
2  DOF/'poi..t  (u^.,  uz) 

CTRIAAX  (triangular  ring) 

generalization  of  TRIARG  which  allows  non-axisymmetric 
deformation  by  expanding  motions  in  Fourier  series 
in  aximuthal  coordinate 


3  DOF/point  ( ,  u,  ,  u  ) 

CTRAPRG  (trapezoidal  ring) 

similar  to  TRIARG  except  trapezoidal  shape 


C3 
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CTRAPAX  (trapezoidal  ring) 

similar  to  TRIAAX  except  trapezoidal  shape 

Miscellaneous  Structural  Elements 
CONMi  (concentrated  mass) 

allows  input  of  a  6x6  symmetric  mass  matrix  at  a 
grid  point 

CDUMi  (dummy  element) 

an  element  for  which  the  user  has  written  his  own 

FORTRAN  subroutines  and  incorporated  them  into  NAS TRAN 

Scalar  Elements  (connect  DOF's  rather  than  grid  points) 

CELASi  (scalar  spring)  •— 

CMASSi  (scalar  mass) 

CDAMPi  (scalar  damper)  • — -j  | — • 

Rigid  Elements 

CRIGD1  (rigid  element) 

all  6  DOF  of  dependent  grid  points  are  coupled  to  all 
6  DOF  of  the  reference  grid  point 
CRIGD2  (rigid  element) 

selected  DOF  of  dependent  grid  points  are  coupled  to 
all  6  DOF  of  the  reference  grid  point 
CRIGD3  (general  rigid  element) 

selected  DOF  of  dependent  grid  points  are  coupled  to 
6  selected  DOF  of  one  or  more  reference  grid  points 
CRIGDR  (rigid  rod) 

a  rod  which  is  rigid  in  extension/compression 

Non-structural  Elements 

CAXIFi  (axisymmetric  fluid  element) 

CFLUIDi  (axisymmetric  fluid  element) 

CSLOTi  (acoustic  cavity  slot  element) 

CHBDY  (heat  transfer  boundary  element) 
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Miscellaneous 

CNGRNT  (identical  elements  indicator) 

designates  secondary  elements  identical 
element  to  avoid  regeneration  of  the 
mass  matrices 


to  a  primary 
stiffness  and 
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TABLE  9-1 

-  Element  Summary 

1-D 

2-D 

3-D 

CBAR 

CTRMEM 

CTETRA 

CROD 

CQDMEM 

CHEXAi 

CON ROD 

CQDMEM1 

C WEDGE 

CTL'BE 

CQDMEM2 

CIHEXi 

CVISC 

CTRBSC 

CTRPLT 

CQDPLT 

CTRIAi 

SCALAR 

CQUADi 

RIGID 

CELASi 

CTRIM6 

CRIGDi 

C MAS Si 

CTRPLT1 

CRIGDR 

CDAMPi 

CTRSHL 

CSHEAR 

CTWIST 

MI  SC.  STRUCT. 

MISC. 

CONMi 

CDUMi 

CNGRNT 

AXI-SYM. 

CCONEAX 

CTORDRG 

CTRIARG 

CTRIAAX 

CTRAPRG 

CTRAPAX 


NON-STRUCT . 
CAXIFi 
CFLUIDi 
CSLOTi 
CHBDY 
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THE  APPLICATION  OF  STRUCTURAL  SYMMETRY 
IN  FINITE  ELEMENT  ANALYSIS 


Summary .  A  brief  review  is  presented  of  the  fundamental  concepts  involved 
in  the  systematic  application  of  structural  symmetry  in  finite  element 
analysis . 


BACKGROUND 


Since  structural  analysts  are  frequently  called  upon  to  analyze 
structures  possessing  symmetry,  it  is  essential  that  the  fundamental 
concepts  of  syimetry  be  sufficiently  well  understood  that  synmetry  can  be 
exploited  systematically  and  with  confidence  [1,2].  The  motivation  for 
wanting  to  exploit  symmetry  is  clear:  when  symmetry  is  present,  the 
engineer  need  model  only  a  portion  of  the  structure,  thereby  saving  both 
his  time  and  the  computer's  time,  with  the  former  probably  being  the  more 
valuable.  For  example,  a  structure  possessing  one  plane  of  symmetry  can 
be  analyzed  by  modeling  only  one-half  of  the  structure,  whether  the  loads 
are  symmetric  or  not.  Even  with  nonsymmetric  loads,  in  which  case  the 
half-structure  would  have  to  be  analyzed  twice,  the  analyst  still  benefits, 
since  a  half-structure  generally  costs  much  less  than  half  as  much  to 
analyze  as  the  complete  structure  would. 


SOME  MOTIVATING  EXAMPLES 


Consider  first  the  following  two-dimensional  example  of  a  simply- 
supported  beam 
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l 
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3  & - *  i 

-a 


Figure  1 
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i  u-.: 


The  symmetry  present  here  is  probably  obvious,  so  that  the  analyst  who 
wanted  to  model  only  one-half  of  the  beam's  span  would  determine  (probably 
by  inspection)  that  the  synmetry  boundary  conditions  to  impose  at  mid-span 
are 

u1  =  0,  R3  =  0  (1) 

(for  the  2-D  problem)  where  u-j  and  are  components  of  the  general  3-D 
displacement  and  rotation  vectors 

u  =  1  +  u2  j  +  u3  k 

(2) 

R  =  8}  i  +  R2  i  +  R3  k 

(It  is  assumed  here  that  grid  points  possess  six  degrees  of  freedom  (DOF). 
Generalization  to  situations  with  mo^e  DOF  per  node  presents  no  problem.) 
Only  slightly  less  obvious  than  the  situation  in  Figure  1  is  the 


Figure  2 

in  which  the  load  is  now  antisymmetric.  In  this  case,  the  anti synme try 
boundary  condition  to  be  applied  at  mid-span  is 

u2  =  0  (3) 

a  condition  which  many  analysts  would  probably  arrive  at  by  inspection. 
Consider  now  the  following  two-dimensional  beam  structure 


Figure  3 
2 
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Like  the  preceding  two  examples,  this  problem  can  also  be  solved  by 
modeling  only  half  the  structure  and  applying  the  appropriate  synmetry 
boundary  conditions.  However,  unlixe  the  preceding  two  examples,  the 
reliance  on  intuition  alone  in  the  application  of  symmetry  would  probably 
fail.  Thus,  what  is  needed  is  a  systematic  procedure  to  follow  with  regard 
to  symmetry.  The  following  sections  sunmarize  such  a  procedure. 

TYPES  OF  SYMMETRY 

In  structural  mechanics,  the  most  corrmonly  encountered  types  of 
symmetry  are  planes  of  symmetry,  axes  of  synmetry,  and  centers  of 
symmetry.  Each  type  of  symmetry  is  characterized  by  some  symmetry 
operation  (reflection,  rotation,  etc.)  which  can  transform  the  structure 
into  an  equivalent  configuration  [2],  The  synmetry  operation  which 
characterizes  a  plane  of  symmetry  is  reflection  in  a  plane,  as,  for  example, 
in  Figure  1.  An  axis  of  symmetry,  exemplified  by  the  structure  in 
Figure  4, 


Figure  4 

is  characterized  by  a  rotation  about  an  axis.  In  this  case,  a  rotation 
of  120°  would  transform  the  structure  and  its  loads  to  an  equivalent 
configuration. 

The  third  type  of  symmetry,  the  center  of  symmetry,  has  as  its 
symmetry  operation  an  inversion  through  a  center  (e.g.,  Figure  5).  In 
Figure  5,  the  center  is  labeled  0. 

We  observe  that  the  above  specifications  of  the  symmetry  operations 
characterizing  a  particular  structure's  symmetry  are  not  necessarily 
unique.  For  example,  we  could  also  characterize  the  symmetry  of 
Figure  5  as  a  sequence  of  two  reflections,  one  in  the  2-3  plane  followed 
by  one  in  the  1-3  plane,  or  vice  versa.  Figure  5  also  serves  as  an 
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Figure  5 

example  of  a  structure  with  an  axis  of  symmetry  (with  a  rotation  angle  of 
180°). 

In  general,  the  plane  of  symmetry  can  be  viewed  as  the  fundamental 
type  of  synrietry,  since  it  can  be  shown  that  all  synmetry  transformations 
of  finite  figures  in  3-0  reduce  to  successive  reflections  in  not  more 
than  three  planes  (which  might  not  even  be  syrmetry  planes) [3]. 

The  identification  of  the  symmetry  possessed  by  some  structure 
requires  not  only  geometrical  symmetry  but  also  symmetry  with  respect  to 
material  properties.  For  example,  if  the  beam  in  Figure  1  were  made  of 
steel  on  the  left  half  and  aluminum  on  the  right,  there  would  be  no 
symmetry  to  exploit.  In  general,  many  other  properties  may  also  play  a 
role  -in  deciding  the  presence  of  symmetry  for  some  problems  (e.g.,  thermal 
radiation  problems  might  require  synmetry  with  regard  to  color). 

Finally,  although  loads  have  been  included  in  some  of  the  examples 
above,  the  application  of  symmetry,  as  will  be  seen,  does  not  depend  on 
the  loads'  being  symmetric.  Thus,  in  determining  the  symmetry  possessed 
by  some  structure,  only  the  structure  (in  the  absence  of  loads)  need  be 
considered.  In  this  case,  the  identification  of  the  symmetry  possessed 
by  a  structure  is  generally  merely  a  matter  of  inspection. 

LOADS 

Once  the  symmetry  properties  of  a  structure  are  identified,  the 
loads  can  be  addressed.  The  question  of  whether  or  not  a  given  system  of 
loads  is  symmetric  or  not  depends  on  the  structure  to  which  it  is  applied. 
Specifically,  a  system  of  loads,  when  applied  to  a  structure  possessing 
certain  symmetry,  is  defined  as  synmetric  if  it  is  brought  into  an 
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equivalent  configuration  by  the  symmetry  operations  of  the  structure. 

The  system  of  loads  is  defined  as  antisymmetric  if  the  symmetry  operations 
plus  a  negation  of  signs  of  all  loads  brings  them  into  an  equivalent 
configuration.  For  example,  the  loads  of  Figures  1,  4,  and  5  are  symmetric 
the  load  system  in  Figure  2  is  anti synmet r i c ,  the  load  system  in  Figure  3 
is  nonsymmetric  (i.e.,  neither  symmetric  nor  ant  i  syrrmetr  i  r ) . 

In  general,  any  nonsymmetric  system  of  loads  can  always  be 
uniquely  decomposed  into  the  sum  of  a  symmetric  and  an  antisymmetric  system 
of  loads  (e.g..  Figure  6). 


Figure  6 


In  Figure  6,  the  symmetric 
are  given  by 


part  of  the  load  F$  and  the  antisymmetric  part 


Fs  *  ?<F1  *  V 
Fa  ’  ?<F,  -  V 


(4) 


where  points  1  and  2  are  image  points  of  each  other. 


THE  GUIDING  PRINCIPLE 


The  principle  upon  which  all  applications  of  symmetry  are  based  is 
that  "equivalent  causes  produce  equivalent  effects,"  or  more  generally, 

"the  effect  is  at  least  as  symmetric  as  the  cause"  [1,4,5].  In  the 
context  of  structural  mechanics,  the  practical  effect  of  this  principle 
is  that  synmetric  loads  produce  symmetric  effects  (displacements,  stresses, 
etc.),  and  antisymmetric  loads  produce  antisymmetric  effects.  In  any 
case,  the  structure  must  be  symmetric  in  order  for  the  principle  to  apply. 
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boundary  conditions 


When  only  a  portion  of  a  symmetric  structure  is  modeled,  the  basic 
principle  provides  the  tool  to  derive  systematical ly  the  symmetry  (or 
antisymmetry)  boundary  conditions. 

Emphasis  here  will  be  restricted  to  planes  of  symmetry  since  the 
symmetry  plane  has  already  been  identified  as  the  fundamental  type  of 
symmetry. 

Consider,  for  example,  the  beam  of  Figures  1  and  2.  We  wish  to 
derive  the  symmetry  and  antisymmetry  boundary  conditions  to  be  applied 
at  the  mid-span  (point  M)  if  only  half  of  the  beam's  span  were  modeled. 

This  is  done  by  (1)  considering  in  turn  each  displacement  component  at  that 
point,  (2)  applying  the  symmetry  (or  antisymmetry)  operations  characterizing 
the  structure  to  that  component  (assumed  to  be  nonzero),  and  (3)  observing 
whether  or  not  that  component  may  in  fact  be  nonzero  and  not  violate 
symmetry.  The  symmetry  operation  for  the  beam  of  Figure  1  is  merely  a 
reflection  into  the  2-3  plane  containing  point  M.  The  antisymmetry 
operations  consist  of  the  same  reflection  followed  by  a  negation  of  sign. 

For  example,  for  the  beam  of  Figure  1,  assume  u^  > 0  at  M.  The  reflection 
results  in  an  image  of  u-|  having  opposite  orientation.  The  additional 
negation  of  sign  (^or  antisymmetry)  yields  a  result  coinciding  with  the 
original  configuration.  Therefore,  u^  must  vanish  at  M  in  order  not  to 
violate  symmetry,  but  u^  may  be  nonzero  for  nonsymmetric  motion. 

Similarly,  we  find  that  u,,  and  u^  vanish  at  M  for  antisymmetry  and  may  be 
nonzero  for  symmetry. 

The  rotational  degrees  of  freedom  require  slightly  different  rules. 
Consider  the  following  problem  which  is  clearly  symmetric.  The  axial 
vector  representation  of  the  symmetric  bending  moments  in  Figure  7  shows 
that  the  reflection  of  an  axial  vector  into  a  plane  requires  an  additional 

(*— i — 0  -  *  i  * 

Figure  7 
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negation  of  sign  compared  to  how  ordinary  vectors  reflect. 

Returning  now  to  the  beam  of  Figure  1,  the  application  of  the 
symmetry  operation  (reflection)  to  the  rotational  components  R^  ,  R^ ,  R^ 
indicates  that,  for  symmetry,  and  R^  must  vanish  in  order  not  to 
violate  symmetry,  and  =  0  for  antisymmetry.  To  summarize,  the  boundary 
conditions  to  impose  at  mid-span  (point  M)  in  Figure  1  are 

u,  =  R0  =  R-.  =  0  for  symmetry 
123  (5) 

R.|  =  u^  =  Uj  =  0  for  antisymmetry 

Since  the  choice  of  coordinate  directions  in  Figures  1  and  2  is 
arbitrary,  the  general i zation  of  the  results  in  eqs.  (5)  is  the  following: 
points  lying  in  a  plane  of  symnetry  can  suffer  no  translation  out  of  the 
plane  and  no  rotation  about  in-plane  lines.  The  anti  symmetry  boundary 
conditions  are  that  the  complementary  degrees  of  freedom  are  constrained. 
The  complementary  nature  of  the  symmetric  and  antisymmetric  boundary 
conditions  is  a  general  result  which  follows  from  the  observation  that  the 
only  distinction  between  anti  syr, me  try  and  symmetry  is  an  additional 
negation  in  the  symmetry  operations. 

NONSYMMETRIC  LOADS 

As  illustrated  in  Figure  6,  any  general  loading  system  can  always 
be  uniquely  decomposed  into  the  sum  of  symmetric  and  antisymmetric 
systems.  For  the  structure  of  Figure  6,  for  example,  the  analyst  would 
model  the  left  half,  say,  and  solve  the  problem  in  two  steps:  (1)  the 
symnetric  part  of  the  load  is  applied  along  with  symmetric  boundary 
conditions  at  M,  and  (?)  the  antisymmetric  part  of  the  load  is  applied 
along  with  antisymmetric  boundary  conditions  imposed  at  M.  Thus,  we  have 
Figure  8,  where  the  symmetric  (S)  and  antisyrmetric  (A)  boundary 


Figure  9 

Taking  the  difference  of  the  S  and  A  solutions  has  the  practical  effect 
of  reversing  the  role  played  by  the  left  and  right  sides.  Thus,  even 
though  only  the  left  side  is  modeled,  the  entire  solution  can  be  obtained. 

MULTIPLE  PLANES  OF  SYMMETRY 


Consider  the  rectangular  region  below  possessing  the  two 
symmetry  planes  indicated: 
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Figure  10 

This  problem  can  be  decomposed  into  four  parts,  as  follows: 

Jl  ,'t  |  t\  U  ,  V  ,  1L  ,1,11 
_ 4 _ —  -H —  +  +  +  ...a;— . 

,  IS  ^  |A  ^  i  S  ^  |A 


i  i 
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Figure  11 
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Here,  the  upper  right  quadrant,  say,  is  modeled,  and  the  four  combinations 
of  symmetry  and  antisymmetry  boundary  conditions  (S-S,  S-A,  A-S,  A-A)  are 
imposed  on  the  points  lying  in  the  two  planes  of  symmetry.  The  four 
solutions  can  be  combined  in  various  ways  to  yield  the  solutions  in  all 
four  quadrants. 


FREE,  UNDAMPED  VIBRATIONS 

The  foregoing  discussion  has  been  devoted  exclusively  to  statics 
problems.  Free,  undamped  vibration  problems  (eigenvalue  problems)  can 
also  exploit  symmetry.  The  calculation  of  all  natural  frequencies  and 
mode  shapes  would  require  one  eigenvalue  run  for  each  unique  combination 
of  symmetry/anti  symmetry  boundary  conditions.  For  example,  the  region  of 
Figure  10,  which  has  two  orthogonal  planes  of  symmetry,  could  be  solved  by 
modeling  only  one  quadrant  and  applying, in  turn,  each  of  the  four 
combinations  of  boundary  conditions. 

It  is  interesting  tj  observe  here  that  the  total  number  of  degrees 
of  freedom  (DOF)  involved  in  the  four  component  problems  of  Figure  10 
exactly  equal  the  original  number  of  DOF  contained  in  the  complete  problem 
[5].  This  follows  as  a  direct  consequence  of  the  symmetry  and  antisynmetry 
boundary  conditions'  involving  complementary  sets  of  DOF.  Thus,  we  have 
"conservation  of  DOF.1  If  this  were  nor  so,  we  would  have  the  disturbing 
situation  in  which  the  mere  application  of  symmetry  results  in  the 
creation  or  destruction  of  DOF .  The  purpose  of  applying  syimetry  is,  of 
course,  to  solve  (with  less  effort)  the  same  problem  rather  than  a 
different  problem. 


TIME-DEPENDENT  PROBLEMS 


All  of  the  preceding  results  for  statics  problems  also  apply  to 
transient  (time-dependent)  situations,  except  that  the  entire  history  of 
time-dependent  loads  must  be  decomposed  into  symmetric  and  antisymmetric 
parts.  This  is  illustrated  in  the  context  of  underwater  shock  response  in 
reference  7. 
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FINITE  ELEMENTS  IN  SYMMETRY  PLANES 

Special  consideration  is  necessary  to  treat  the  situation  in  which 
elements  lie  entirely  in  a  plane  of  symmetry  (i.e.,  the  grid  points  which 
define  the  element  lie  entirely  in  the  plane).  For  example,  in  the 
stiffened  plate  shown  below. 


Figure  12 


the  beam  stiffener  (modeled  with  beam  elements)  lies  entirely  in  the  xz- 
plane,  which  is  a  plane  of  structural  symmetry.  Although  the  symmetry 
boundary  conditions  are  unaffected  by  this  situation,  care  must  be 
exercised  in  computing  the  geometrical  properties  of  a  beam  element  lying 
in  the  symmetry  plane.  In  particular,  the  properties  for  each  "half¬ 
element"  should  be  specified  so  that  the  "half-element"  receives  one-half 
the  total  stiffness  rather  than  one-half  the  cross-section.  For  example, 
properties  such  as  area  A,  cross-sectional  moments  of  inertia  1^,  Ip,  and 
I-jp,  and  torsional  constant  J  would  first  be  computed  for  the  full  cross- 
section  before  entering  one-half  of  those  values.  (Note  that  J  and  one  of 
I-j  or  Ip  do  not  depend  linearly  on  individual  cross-sectional  dimensions). 

To  prove  the  validity  of  the  above  approach,  we  need  only  treat 
each  half  of  the  symmetric  structure  as  a  "super  element"  involving  many 
grid  points.  Then,  if  the  two  sides  were  to  be  recombined  using  the  usual 
rules  of  matrix  assembly,  the  resulting  stiffness  matrix  would  have  to  be 
the  correct  stiffness  matrix  for  the  entire  structure.  Thus,  when  an 
element  is  cut  in  half  by  a  symmetry  plane,  each  side  receives  one-half 
of  the  total  stiffness. 
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A  RETURN  TO  FIGURE  3 


The  third  motivating  example  used  earlier  (Figure  3)  is  now  seen 
as  a  routine  application  of  symnetry.  Let  us  characterize  the  symmetry  of 
the  structure  as  the  sequence  of  two  reflections,  one  in  the  1-3  plane 
containing  the  center  0  and  one  in  the  2-3  plane  containing  the  center. 

By  considering  in  turn  each  of,  say,  six  DOF  at  the  center  point,  we  find 
that  the  displacement  vector  at  point  0  must  satisfy 


u-j  =  u^  =  R^  =  R2  =  0  for  symmetry 

u^  =  R^  =  0  for  antisymmetry 


Thus  the  problem  can  be  decomposed  as  follows: 


^4* 


(6) 


h 


Figure  13 

where  the  S  and  A  boundary  conditions  are  given  in  eqs.  (6). 
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GRID  POINT  SEQUENCING  CONSIDERATIONS 
IN  FINITE  ELEMENT  ANALYSIS 


Summary.  A  brief  review  of  the  definitions  of  matrix  bandwidth,  profile, 
and  wavefront,  and  their  implications  in  finite  element  analysis,  is 
presented. 


BACKGROUND 

Many  problems  of  scientific  and  engineering  interest  reduce  to  the 
numerical  problem  of  having  to  solve  a  large  set  of  linear  algebraic 
equations  such  as  (in  matrix  form) 

Ax  =  b  (1) 

where  the  vector  b  and  the  square  matrix  A  are  known,  and  the  unknown 
vector  x  is  sought.  In  finite  element  and  other  applications,  A  is  also 
sparsely  populated  (i.e.,  it  contains  far  more  zeros  than  nonzeros), 
since  the  procedure  under  which  finite  element  matrices  are  assembled 
dictates  that  the  off-diagonal  matrix  terms  coupling  any  two  degrees  of 
freedom  to  each  other  are  zero  unless  those  degrees  of  freedom  are 
common  to  the  same  finite  element.  It  also  follows  that  the  locations 
of  the  nonzero  matrix  elements  of  the  matrix  A  depend  solely  on  the 
ordering  of  the  unknowns.  Thus,  it  is  possible  with  sparse  matrices  to 
choose  an  ordering  which  results  in  the  nonzeros'  being  located  in  such  a 
way  as  to  be  convenient  for  subsequent  matrix  operations  such  as  equation 
solving  or  eigenvalue  extraction.  Many  such  algorithms  have  been 
expressly  written  to  operate  very  efficiently  on  matrices  possessing  small 
bandwidth,  profile,  or  wavefront.  This  is  accomplished  by  avoiding 
arithmetic  operations  on  matrix  elements  known  to  be  zero.  As  a  result, 
the  execution  time  for  a  "band  solver",  for  example,  would  be  0(NB:)  for 
large  N  and  B, where  N  is  the  matrix  order  and  B  the  bandwidth.  For  a 
given  structural  model,  N  is  fixed,  but  B  depends  on  the  ordering  of  the 
unknowns  (grid  points).  Clearly,  in  this  case,  it  is  desirable  to  reduce 
B  as  much  as  possible. 
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DEFINITIONS 


Although  the  definitions  to  be  given  here  are  reasonably  standard 
(at  least  in  finite  element  circles),  uniformity  of  definitions  and 
notation  among  the  various  workers  in  the  field  does  not  yet  exist. 

Given  a  symmetric  square  matrix  A  of  order  N,  we  define  a  "row 
bandwidth"  b-  for  row  i  to  be  the  number  of  columns  from  the  first  nonzero 
in  the  row  to  the  diagonal,  inclusive.  Numerically,  b^  exceeds  by  unity 
the  difference  between  i  and  the  column  index  of  the  first  nonzero  entry  of 
row  i  of  A.  Then  the  matrix  bandwidth  B  and  profile  P  are  defined  as 


max  b . 
ifN  1 


(2) 


P  = 


N 

I 

i  =  l 


b. 

■ i 


(3) 


Let  c.j  denote  the  number  of  active  columns  in  row  i.  By  definition, 
a  column  j  is  active  in  row  i  if  i  2  i  and  there  is  a  nonzero  entry  in 
that  column  in  any  row  with  index  k  S  i.  The  matrix  wavefront  W  is  then 
defined  as 


iSN  ' 

Sometimes  c-  is  referred  to  as  the  row  wavefront  for  row  i.  Since  the 
matrix  A  is  symmetric, 

N  N 

P  =  I  b.  =  I  c.  (5) 

i  =  l  1  1=1  1 


The  wavefront  W  is  sometimes  called  the  maximum  wavefront  W _  to 

J  nid  X 

distinguish  it  from  the  average  wavefront  snd  root -mean-square 

wavefront  W  defined  as 
— — — - —  rms 

1  N  P 

Wavg  =  N  iJ1  Ci  =  N 


(6) 


W 


rms 


(7) 
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As  a  consequence  of  the  above  definitions,  it  follows  that,  for  a  given 
matrix , 


U  ^  Uc  -BIN 

avg  rms  max 


(8) 


The  first  two  inequalities  would  be  equalities  only  for  uninteresting 
special  cases  such  as  diagonal  matrices. 

We  define  the  degree  d^  of  node  i  to  be  the  number  of  other  nodes 
to  which  it  is  connected;  i.e.,  to  be  more  precise,  is  the  number  of 
nonzero  off-diagonal  terms  in  row  i  of  the  matrix  A.  (This  implies,  for 
example,  that  diagonally  opposite  nodes  in  a  quadrilateral  element  are 
"connected"  to  each  other.)  Hence,  the  maximum  nodal  degree  M  is 

M  =  “jj  d.  (9) 

The  number  of  unique  edges  E  is  defined  to  be  the  number  of  nonzero  off- 
diagonal  terms  above  the  diagonal.  Hence,  for  a  symmetric  matrix. 


1  N 

E  =  j  I  d.  (10) 

c  i  =  1  1 

Thus  the  total  number  of  norizeros  in  A  is  2E+N,  and  the  density  o  of  the 
matrix  A  is 


N 

Note  that,  in  the  above  definitions,  the  diagonal  entries  of  the 
matrix  A  are  included  in  b-  and  c.  (and  hence  in  B,  P,  W  ,  W  ,  and 

1  1  IDa  X  a  VQ 

W  ).  Therefore,  these  definitions  differ  from  those  of  some  authors, 
rms 

but  conform  to  those,  for  example,  in  the  NASTRAN  literature.  Except  for 

the  rms  wavefront  W  ,  it  is  easy  to  convert  the  various  parameters  from 
rms 

one  convention  (including  the  diagonal)  to  the  other  (not  including  the 
diagonal ) . 

Also  note  that,  in  this  context,  the  order  N  of  the  matrix  A  is 

sometimes  taken  to  be  the  same  as  the  number  of  nodes.  In  general  finite 

element  usage,  however,  each  node  (grid  point)  has  several  degrees  of 

freedom  (DOF),  not  just  one.  For  structures  having,  say,  six  DOF  per 

node,  the  actual  DOF  values  of  B,  W  ,  W  ,  or  W _  would  be  (in  the 

max  avg  rms 

absence  of  constraints)  six  times  their  corresponding  grid  point  values. 
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Example 

The  above  definitions  ( 2 ) - ( 7 )  can  be  illustrated  by  the  following 
simple  example.  In  Figure  1  is  shown  a  matrix  of  order  six.  In  each  row 
and  column  a  line  is  drawn  from  the  first  nonzero  to  the  diagonal.  Thus 
b-  is  the  number  of  columns  traversed  by  the  solid  line  in  row  i. 

-X  X  ’  -:i 

X  i  X  < 


bj. 

1 
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3 
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X-  — 


i 


£i 

3 

5 

4 
3 
2 
1 


£l 

9 

25 

16 
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4 
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E  =18  E  =  64 


Figure  1  -  Example  illustrating  definitions  of  matrix 
bandwidth,  profile,  and  wavefront. 

Similarly,  the  number  of  active  columns  c.  in  row  i  is  the  number  of 

vertical  lines  in  row  i  to  the  right  of  and  including  the  diagonal. 

Thus,  from  the  above  definitions,  B=6,  W  =5,  P=18,  W  =3.0,  and  W  =3.3. 
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THE  RELATIONSHIP  TO  STRUCTURES 


Consider  the  one-dimensional  six  DOF  system  of  six  scalar  springs 


shown  below. 


£  J 

] — vVW — t — v/WV — ® — AAV - 14  5  6 

v _ — _  _ j — AM/  1  ° — v'VW'  0 


-\MV 


Figure  2 

For  each  spring,  the  element  stiffness  matrix  is 

’  1  -1 
-1  1 


(12) 


where  k  is  the  spring  stiffness.  For  the  nodal  numbering  indicated  in 


1 1  -  > 


Figure  2,  the  6x6  system  stiffness  matrix  looks  like 


where  an  X  indicates  the  location  of  a  nonzero  entry.  From  equation  (13) 
and  the  definition  (2),  the  matrix  bandwidth  B  is  4.  Observe  that  the 
bandwidth  can  also  be  obtained  directly  from  the  structure  (Figure  2)  by 
adding  unity  to  the  maximum  numerical  difference  between  connected  node 
numbers  (node  1  is  connected  to  node  4).  The  same  result  is  also 
obtained  for  the  following  numbering: 


Figure  3 


This  is  because,  from  the  point  of  view  of  the  matrix  connectivity,  there 
is  no  difference  between  the  structures  in  Figures  2  and  3,  since  the 
ordering  of  the  unknowns  is  the  same.  Some  structural  programs  (e.g.  , 
NASTRAN)  allow  the  user  to  specify  grid  point  numbers  as  in  Figure  3, 
rather  than  consecutively  from  1  to  N.  However,  in  order  to  compute  the 
matrix  bandwidth  directly  from  this  structure  (by  looking  at  the  maximum 
numerical  difference  between  connected  node  numbers).  Figure  3  would 
first  have  to  be  simplified  to  Figure  2. 

To  illustrate  the  difference  that  sequencing  makes,  consider  instead 
the  following  numbering 


Figure  4 


5 


1 1  -i. 


Here  the  bandwidth  is  6,  so  that  the  ordering  of  Figure  2  is  to  be 
preferred  over  that  of  Figure  4.  However,  a  still  better  sequence  (i.e., 
one  with  a  smaller  bandwidth)  is 


where  B=3. 

The  same  concepts  can  also  be  applied  to  two-  and  three-dimensional 
structures.  For  example,  consider  the  plate  below  modeled  with  a  2x4 
array  of  quadrilateral  elements. 
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Figure  6 

Here  the  grid  point  bandwidth  is  7  (recall  that  all  nodes  in  a  given 
element  are  "connected"  to  all  other  nodes  in  the  same  element).  A  better 
sequence  (i.e.,  one  with  a  smaller  bandwidth)  would  number  first  across 
the  "short"  direction  ("short"  in  the  sense  of  number  of  nodes  rather  than 
actual  .distance) : 
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Figure  7 

With  this  sequence  the  grid  point  bandwidth  is  now  5. 

In  general,  the  plates  of  Figures  6  or  7  have  more  than  one  DOF  per 
node.  Thus,  although  it  suffices  to  consider  only  grid  point  bandwidth 


6 


II-/ 


when  picking  an  ordering,  the  actual  DOF  bandwidth  which  the  equation 
solver  encounters  would  be  much  larger.  For  example,  structures  having 
6  DOF/node  and  a  grid  point  bandwidth  of  B  would  have  a  DOF  bandwidth  of 
6B . 

Although  the  above  discussion  was  written  from  the  point  of  view  of 
matrix  bandwidth,  similar  comments  could  be  made  instead  from  the  point 
of  view  of  matrix  profile  or  wavefront.  For  e*ample,  NASTRAN's  level  16 
contains  a  decomposition  routine  which  operates  fastest  on  those  matrices 
having  smallest  rms  wavefront. 


AUTOMATIC  RESEQUENCERS 


Although  the  preceding  sections  define  the  various  terms  and  show 
how  one  might  compute  the  bandwidth,  profile,  or  wavefront  for  a  given 
matrix,  such  calculations  would  clearly  be  very  tedious  for  all  but  the 
smallest  structures.  Even  more  difficult,  in  general,  is  the  job  of 
resequencing  the  nodal  labels  to  reduce  the  parameter  of  interest.  This 
is  especially  true  for  large,  complicated  meshes  or  those  generated 
automatically  on  a  computer. 

Fortunately,  a  large  number  of  algorithms  have  been  developed  to 
automate  the  assignment  of  grid  point  labels,  given  the  connectivity  of 
the  mesh.  Since  it  is  clearly  impractical  to  check  each  of  the  N! 
possible  sequences  associated  with  a  given  matrix  A  of  order  N,  each 
algorithm  attempts  some  (presumably)  rational  strategy  for  arriving 
quickly  at  a  grid  point  sequence. 

For  NASTRAN,  for  example,  two  preprocessors  are  currently  being  used 
to  resequence  nodes:  BANDIT  [1-3],  which  uses  the  Cuthill-McKee  [4]  and 
Gi bbs-Pool e-Stockmeyer  [5]  algorithms,  and  WAVEFRONT  [6],  which  uses  the 
Levy  algorithm  [7,8].  Some  good  comparisons  of  several  resequencing 
algorithms  have  been  made  by  Cuthill  [9]  and  Gibbs,  et  a]  [10]. 
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1.  BANDIT  version  :  ve r-. i on  9,  updated  ; .  i  "v 

2.  Compatible  NASTRAN  Versions:  NASA’ s  level  IT  and  below),  MSt  NASTRAN, 

Navv  NASTRAN 

3.  I nput : 

a.  max invcn :  a  standard  NASTRAN  data  deck  ;  ID  or  NASTRAN  thru  KM'DATA) 

plus  $  option  cards,  if  any 

b.  minimum:  $  option  cards,  BEGIN  BULK,  element  connection  cards, 

LNT'DA'.  A 

4.  Output : 

a.  printed  output 

b.  punched  output  (SEQGP  cards  or  entire  deck) 

c.  file  (unit  81  containing  complete  deck  plus  SEQGP  cards,  this  file 

is  suitable  to  be  used  as  input  to  NASTRAN 

5 .  Elements  Recognized: 


CELAS1 

CELAS2 

CDArtPl 

CDAfIPe 

CHASS1 

CPIASS2 

CROD 

CTUBE 

CWISC 

CDANP3 

CDAHP4 

CELAS3 

CEIAS4 

CWASS3 

CHASS4 

CAXIF2 

CAXIF3 

CAXIF4 

CBAR 

CCONEAX 

CFLUID2 

CFIUID3 

CFLUID4 

CHBDV 

CHEXA1 

CHEXA2 

CHTTRI2 

CIS2D4 

CIS2D8 

CIS3D8 

CIS3D2® 

C0NPI1 

C0NH2 

CONROD 

CODNEN 

CQDMEP11 

CQDHEN2 

CQDPLT 

CQUADl 

CQUAD2 

CSHEAR 

CSL0T3 

CS10T4 

ctetra 

CTORDRG 

CTRAPRG 

CTRBSC 

CTRIA1 

CTRIA2 

CTRIARG 

ctruen 

CTRPLT 

CTUIST 

CUEDGE 

CDunnv 

cduhi 

CDUH2 

CDUPI3 

CDUN4 

CDimS 

CDUH6 

CDUH7 

CDUH8 

CDUH9 

CTRIAX6 

CTRIH6 

CDANP4* 

CELAS4* 

CHASS4* 

CDAHP2* 

CELAS2* 

CNASS2* 

CONFU* 

C0NH2* 

CONROD* 

CIHEX1 

CIHEX2 

CIHEX3 

CTRAPAX 

CTRIAAX 

CQUADTS 

CTRIATS 

CQDflEW3 

CHEX8 

CMEX20 

CTRPLT 1 

CTRSHL 

CR1GD1 

CRIGD2 

CRIGDR 

CBEAN 

CTPIA3 

CFTUBE 

CHEXA 

cpehta 

CQUAD4 

TM- 184  "-03 


6.  Reduction  Approach:  Uses  Cuthill-McKee  (CM)  and/or  Gibbs- Poole -Stockmeyer 

(GPS)  methods  to  reduce  matrix  bandwidth,  profile, 
wavefront,  or  rms  wavefront. 


7 .  Core 


ruirements : 


Total  Core  =  Program  +  Working  Storage 

where  Program  =  4 7 Kg  words  on  CDC 
=  14 5K  bytes  on  IBM 
=  24K  words  on  UNIVAC,  Honeywell 

Working  Storage  Required  =  +  8)N 

where  N  =  number  of  grid  points 

M  =  maximum  nodal  degree  (the  maximum  number  of  nodes  connected 
to  any  node) 

MV  =  integer  packing  density  (integers/word) 

Working  Storage  =  open  core  on  open  core  versions  of  BANDIT 

CDC:  MV  =  6  for  N  s  510 

5  for  510  <  N  -  2045 
4  for  2045  <  N  <  16380 
3  for  16380  <  N  i  524286 


IBM:  MV 


N  1  32766 


'rNTYAC ,  Honeywell:  NW  =  4  for  N  I  508 

3  for  508  <  N  4095 


2  for  4095  <  N  5  262142 


)tion  Cards: 


location:  anywhere  before  BEGIN  BULK 
general  format:  $ KEYWORD 1  KEYWORD 2 


c.  rules: 


(1)  S  in  column  1 

(2)  KEYW0RD1  starts  in  column  2 

(3)  keywords  separated  by  one  or  more  blanks 

(4)  no  embedded  blanks  in  keywords 

(5)  the  first  two  letters  of  each  keyword  are  required  for 
recognition 
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9.  $  Option  Cards:  (default  underlmedj 

a.  For  General  Use: 

$ SEQUENCE  (NO,  YES)  Is  resequencing  to  be  performed? 

$GRID  N  Upper  bound  on  number  of  grid  points. 

This  card  must  be  used  with  solid  elements  and  MFC's  since 
default  M  (maximum  nodal  degree)  is  about  19.  It  is 
recommended  to  use  it  with  all  runs  since  it  is  used  for 
efficient  allocation  of  core. 

$ CONFIG  N  Computer  model  (from  NASI RAN  manual). 

Used  in  estimating  NASTRAN  decomposition  time. 

SCONTIG  N,M,L  N  =  computer  model  (from  NASTRAN  manual) 

M  =  computer  for  which  decomposition  time 
estimate  is  desired  if  different  from 
one  RANDI'l  is  on  (M=l  for  CDC, 

2  for  IBM,  3  for  UNTYAC). 

I.  =  flag  to  request  printout  of  all 
NASTRAN  multiply-add  time  constants 
(0  =  no,  1  =  yes) 

$ PUNCH  (NONE ,  SFQGP ,  ALL)  What  should  be  punched? 

SCRITERION  (BAND,  PROFILE,  WA\T  FRONT ,  RMS)  What  should  be  reduced1? 
Recommendations : 

BAND  for  NASTRAN  Level  15.5  and  below 

RMS  for  NASTRAN  Level  15.9  and  above  and  MSC  NASTRAN 

$MFTH0D  (CM,  GPS,  BOTH)  By  what  method? 

SMPC  (NO,  YES)  Take  MPC's  into  account? 

"YES"  generates,  for  each  MR?  equation  in  deck,  additional 
connections  between  the  independent  points  and  every  other  point 
to  which  the  dependent  point  is  connected.  Dependent  points  can 
be  eliminated  from  connection  table  by  using  SIGNORE. 

SPRINT  (MIN,  MAX)  What  printed  output? 

"MIN"  is  adequate  for  most  purposes.  "MAX"  generates 
additional  connection  tables  and  nodal  lists. 

SIGNORI;  G1,G2,...  Grid  points  to  ignore. 

Nodes  ignored  are  eliminated  from  the  connection  table  and 
sequenced  last.  This  should  be  used,  for  example,  for  nodes 
of  very  high  degree  compared  to  other  nodes  in  the  structure. 

$ADD  N  Add  N  to  new  sequence  numbers. 

May  be  used  to  avoid  duplication  of  internal  numbers  if  not 
all  nodes  of  a  structure  are  being  sequenced  in  one  run. 
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$ELtMENTS  (NO,  YES)  List  BANDIT'S  element  library? 

$APPEND  CNAME  NCON  IFLD  User-defined  connection  card. 

CNAME  =  name  of  connection  card  (e.g.,  CBAR)  left-adjusted 
starting  in  column  9. 

NCON  =  number  of  connections  on  card  (i.e.,  nodes  in 
element)  ?  1. 

IFLD  =  NASTRAN  field  number  on  parent  card  in  which  first 
connection  appears  s  9. 

NCON  and  IFLD  may  appear  anywhere  in  columns  17-32  separated 
by  one  or  more  blanks.  No  long- field  connection  cards  may 
be  defined.  Connections  most  be  listed  consecutively  on 
parent  and  continuation  cards,  if  any.  Each  SAPPEND  card 
defines  a  new  element  type. 


b. 


For  Particular  Users: 

$ NASTRAN  (NO,  YES)  NASTRAN  to  follow  BANDIT0  (IBM) 

"YES"  generates  a  condition  code  5  after  a  successful  completion. 
$  INSERT  Location  of  cards  to  insert. 

$  INSERT  N  Number  and  location  of  cards  to  insert. 


May  be  used  by  remote  users  to  insert  checkpoint  dictionary’ 
from  disk  file  into  executive  control  deck. 


$ LINES  N 
$PLUS  ♦ 

Allows  user  to 
$ DIMENSION  N 
SHI CORE  N 


Number  of  printed  lines  per  page. 
User-defined  plus  sign. 

input  his  own  special  plus  sign,  if  necessary. 
Dimension  of  a  scratch  array. 

Amount  of  core  requested  in  words.  (UNI VAC) 


c. 


For  Program  Developer  : 
STABLE  (NO,  YES) 

SSTART  G1 ,GZ , . . . 

$ DEGREE  N 
SSPRING  (NO,  YES) 


Output  connection  table? 
User-supplied  CM  starting  nodes. 
Ignore  nodes  of  degree  exceeding  N. 
Generate  scalar  springs? 


The  springs  (CELAS3)  have  same  connectivity  as  original 
structure. 


10.  Installation-dependent  Remarks: 

a.  On  CDC  machines,  the  automatic  reduction  of  field  length  at 
execution  time  should  be  suppressed,  e.g.,  with  an  RFL  card. 

b.  Unless  modified  locally,  IBM  and  Honeywell  versions  are  not 
open  core  programs,  but  fixed  core.  Hence,  calls  by  BANDIT  for 
more  core  require  the  change  of  two  statements  in  the  main  program. 
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11.  Additions  to  BANDIT  version  9: 


a.  Open  core  and  HICORL  on  UN I  VAC . 

b.  Eliminate  backspace  of  Unit  5  on  IBM  and  Honeywell. 

c.  Min.  nodal  degree  printout  in  summary. 

d.  User-selected  CM  starting  nodes  fix. 

e.  Case  Control  card  counter. 

f.  New  Level  17  configurations  and  time  constants. 

g.  Subroutine  READIT  efficiencies. 

h.  Recovery'  of  SEQGP  cards  generated  by  CM  if  abort  in  GPS  due 
to  exceeding  scratch  dimension. 

i.  $APPEND  card  to  define  connection  card  at  execution  time. 

j.  CRIGD1  with  THRU  option. 

k.  Warning  for  illegal  ENDDATA  format. 

l.  Optional  printout  of  multiply-add  time  constants. 

m.  Reset  of  SDIMENSION  value  if  $GRID  declared. 

n.  Time  and  disk  space  efficiencies  with  MPC  equations. 
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DTNSRDC  ISSUES  THREE  TYPES  OF  REPORTS 

1.  DTNSRDC  REPORTS,  A  FORMAL  SERIES,  CONTAIN  INFORMATION  OF  PERMANENT  TECH 
NICAL  VALUE.  THEY  CARRY  A  CONSECUTIVE  NUMERICAL  IDENTIFICATION  REGARDLESS  OF 
THEIR  CLASSIFICATION  OR  THE  ORIGINATING  DEPARTMENT. 

2.  DEPARTMENTAL  REPORTS,  A  SEMIFORMAL  SERIES.  CONTAIN  INFORMATION  OF  A  PRELIM 
INARY,  TEMPORARY,  OR  PROPRIETARY  NATURE  OR  OF  LIMITED  INTEREST  OR  SIGNIFICANCE 
THEY  CARRY  A  DEPARTMENTAL  ALPHANUMERICAL  IDENTIFICATION. 

3.  TECHNICAL  MEMORANDA,  AN  INFORMAL  SERIES.  CONTAIN  TECHNICAL  DOCUMENTATION 
OF  LIMITED  USE  AND  IN  I  EREST.  THEY  ARE  PRIMARILY  WORKING  PAPERS  INTENDED  FOR  IN 
TERNAL  USE.  THEY  CARRY  AN  IDENTIFYING  NUMBER  WHICH  INDICATES  THEIR  TYPE  AND  THE 
NUMERICAL  CODE  OF  THE  ORIGINATING  DEPARTMENT.  ANY  DISTRIBUTION  OUTSIDE  DTNSRDC 
MUST  BE  APPROVED  BY  THE  HEAD  OF  THE  ORIGINATING  DEPARTMENT  ON  A  CASE  BY  CASE 
BASIS. 


